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ABSTRACT 

We present a 3D simulation of the dynamic emergence of a twisted magnetic 
flux tube from the top layer of the solar convection zone into the solar atmosphere 
and corona. It is found that after a brief initial stage of flux emergence during 
which the two polarities of the bipolar region become separated and the tubes 
intersecting the photosphere become vertical, significant rotational motion sets 
in within each polarity. The rotational motions of the two polarities are found 
to twist up the inner field lines of the emerged fields such that they change 
their orientation into an inverse configuration (i.e. pointing from the negative 
polarity to the positive polarity over the neutral line). As a result, a flux rope 
with sigmoid-shaped, dipped core fields form in the corona, and the center of 
the flux rope rises in the corona with increasing velocity as the twisting of the 
flux rope footpoints continues. The rotational motion in the two polarities is a 
result of propagation of non-linear torsional Alfven waves along the flux tube, 
which transports significant twist from the tube's interior portion towards its 
expanded coronal portion. This is a basic process whereby twisted flux ropes are 
developed in the corona with increasing twist and magnetic energy, leading up 
to solar eruptions. 

Subject headings: MHD — Sun: active regions — Sun: coronal mass ejections 
— Sun: flux emergence — Sun: magnetic fields 

1. Introduction 

How twisted magnetic fields emerge from the dense, convectively unstable solar convec- 
tion zone into the stably stratified, rarefied solar atmosphere and corona is a fundamentally 
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important question for understanding the formation of solar active regions and the develop- 
ment of precursor structures for solar eruptions such as flares and CMEs. Pioneered by the 
earlier work of Shibata and collaborators ( Shibata et al.|1989 ), a large body of 3D MHD sim- 



ulations of the emergence of twisted magnetic flux tubes into the solar atmosphere (e.g. Fan 
200T| |Magara and Longcope 1[200T| [20031 |Magara |[2004| [Manchester et al. ||2004t |Archontis 



et al. 2004 Murray et al. 2006 Magara 2006 Archontis and Toeroek 2008) have been 



carried out in recent years, which have provided important insight into the above dynamic 
process. It has been shown that magnetic flux reaching the photosphere can undergo a dy- 
namic expansion into the atmosphere as a result of the non-linear growth of the magnetic 
buoyancy instability. The conditions required of the emergent tubes for the onset of such 



dynamic expansion into the atmosphere has been examined by Archontis et al. (2004) and 



Murray et al. (2006). 



These simulations have also shown that it is difficult for a twisted flux tube to rise 
bodily into the corona as a whole due to the heavy plasma that is trapped at the bottom 



concave (or U-shaped) portions of the winding field lines (e.g. Fan 2001 Manchester et al. 



2004 Archontis et al. 2004; Magara 2006). While the upper parts of the helical field lines 
of the twisted flux tube expand into the atmosphere and the corona, the U-shaped parts 
of the winding field lines remain largely trapped at and below the photosphere layer, and 
the center of the original tube axis ceases to rise a couple of pressure scale heights above 
the photosphere. Despite this, a twisted coronal flux rope with sigmoid-shaped concave 
upturning field lines threading under a central axis that rises into coronal heights is found 



to develop in the end in many simulations (e.g. 


Magara and Longcope 


2001. 


2003 


Magara 


2004 


Manchester et al. 


2004; 


Magara 2006 


Archontis and Toeroek 


2008 


)• 


Manchester 



et al. (2004) found that shear flows on the photosphere driven by the magnetic tension 
force continually transport axial flux and magnetic energy to the upper expanded portion of 
the emerged field, causing it to eventually "pinch off" from the lower mass-laden part via 
magnetic reconnection in a current sheet forming in the lower atmosphere. This reconnection 
produces a coronal flux rope with a new central axial field line that rises into the corona 
with increasing speed. Such development of a coronal flux rope is found to take place more 
readily with a more sharply arched subsurface tube for which the emerging segment breaking 



through the photosphere contains fewer turns of field line twist. Magara and Longcope 



(2003) shows that during the emergence of a twisted flux tube, both the vertical emergence 
of magnetic flux and the horizontal motions at the footpoints of the emerged field lines inject 
magnetic energy and helicity (or twist) into the atmosphere. The contribution by the direct 
emergence dominates at the early phase of flux emergence, while the horizontal flows become 



a dominant contributor of helicity injection later. Magara (2006) studied in more detail the 



horizontal flows on the photosphere during flux emergence and found that as soon as the flux 
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tubes intersecting the photosphere become vertical and the two polarity flux concentrations 
on the photosphere become separated, a rotational flow appears in each of the polarity flux 
concentrations. 

In the present paper, we study further the role of the shear flows and the rotational flows 
in driving the formation of a coronal flux rope during flux emergence. We present a simulation 
of the emergence of a twisted fi-shaped subsurface flux tube whose emerging segment initially 
breaking through the photosphere contains less than one full wind of field-line twist. As in 



the simulation of Magara (2006), we find a brief initial phase of flux emergence during 
which the two polarity flux concentrations on the photosphere undergo a shearing motion 
along the neutral line and become separated. Then a prominent rotational/vortical motion 
develops within each polarity, reminiscent of sunspot rotations that have been observed 



occasionally (e.g. Brown et al. 2003). The rotational motions of the two polarities are due 
to torsional Alfven waves propagating along the twisted flux tube, transporting magnetic 
twist from the tube's interior portion, where the rate of twist is high, to its greatly stretched 
coronal portion, where the rate of twist is low, a result that was predicted by an earlier 



idealized analytical model of Longcope and Welsch (2000). The rotational motions of the 



two polarities are found to persist throughout the later phase of the evolution and steadily 
inject magnetic helicity into the atmosphere. The combined effect of the initial shear flow 
and the subsequent rotational motions of the two polarities is to twist up the inner field 
lines of the emerged field, causing them to rotate from their initial normal configuration (i.e. 
arching over the polarity inversion line from the positive to the negative polarity), into an 
inverse configuration (directed from the negative polarity to the positive polarity over the 
neutral line). In this manner, a flux rope with sigmoid-shaped dipped core fields forms in the 
atmosphere, with the center of the flux rope, as represented by the O-point of the transverse 
field in the central cross-section, rising into coronal heights. This field reconfiguration is a 
major new finding of our study, namely that the horizontal rotation of the field direction of 
the lower emerged flux into an inverse configuration causes an upward propagation of the 
O-point, and initiates the development of the coronal flux rope. Similar to many previous 



Magara and Longcope |2003 


Magara |2004 


Manchester et al. |2004 


Magara 



2006), a current sheet forms in the lower atmosphere, and the field lines whose lower parts 
going through the current sheet all display a sigmoid morphology, which may correspond to 
the observed X-ray sigmoid loops in active regions. 

The remainder of the paper is structured as follows. In §2 we describe the numerical 
MHD model and the detailed setup of the simulation. In §3, we describe the simulation 
results, with focuses on examining the nature of the rotational motion of the two polarities 
on the photosphere and how it drives the formation of a flux rope in the corona. Finally in 
§4 we summarize the conclusions and discuss the implications of our simulation results. 
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The numerical model 



In this study we solve numerically the following ideal magnetohydrodynamic (MHD) 
equations in a Cartesian domain: 



dp 
dt 



V • (pv) = 0, 



dt 
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In the above v, B, p, p, T, R, p, e, e denote respectively the velocity field, the magnetic field, 
density, pressure, temperature, the gas constant, the mean molecular weight, the internal 
energy per unit volume, and the total energy (internal+kinetic+magnetic) per unit volume. 
The gravitational acceleration g = —gz is a constant and is directed in the —z direction. The 
above ideal MHD equations are solved numerically without any explicit viscosity, magnetic 
diffusion, and non-adiabatic effects, although numerical dissipations are present. Thus we 
do not physically treat the radiative energy transfer, thermal conduction, coronal heating, 
and the effects of partial ionization, which all play important roles in the thermal energy 
evolution of the plasma in the solar atmosphere. The thermal properties of the plasma in our 
simulations are therefore unrealistic. Our study here focuses on the dynamic and topological 
properties of the 3D emerging magnetic field, given the above simplifications. 

The basic numerical schemes we use to solve the above MHD equations are as follows. 



The equations are discretized in space using a staggered finite-difference scheme (Stone 



and Norman 1992a), and advanced in time with an explicit, second order accurate, two- 
step predictor-corrector time stepping. A modified, second order accurate Lax-Friedrichs 



scheme similar to that described in Rempel, Schussler, and Knolker (2008, see eq. (A3) in 



that paper) is applied for evaluating the fluxes in the continuity, momentum, and energy 
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equations. Compared to the standard second order Lax-Friedrichs scheme, this scheme 
significantly reduces numerical diffusivity for regions of smooth variation, while retaining 
the same robustness in regions of shocks. A method of characteristics that is upwind in 



the Alfven waves (Stone and Norman 1992b) is used for evaluating the v x B term in the 



induction equation, and the constrained transport scheme is used to ensure V • B = to 
machine round-off errors. 

The equations are non-dimensionalized using the following normalizing constants: p p h ~ 
2.7 x 10~ 7 g cm -3 (photospheric density), if p h ~ 150 km (pressure scale height at the pho- 
tosphere), T p h « 5100 K (photospheric temperature), v t h = (7£T p h/ /i) 1 / 2 ~ 6.5 km s -1 
(isothermal sound speed at the photosphere), t ph = if ph /t> th « 23 s, p p h = Pph^th (pho- 
tospheric pressure) and B ph = (Airpp^v^) 1 ^ 2 ~ 1200 G, as the units for density, length, 
temperature, velocity, time, pressure, and magnetic field respectively. For the remainder 
of the paper, quantities and equations are expressed in the above units unless otherwise 
specified. 

The setup of the numerical simulation is illustrated in Figure 1. The simulation domain 
(top panel of Figure 1) is (—69.3,-63.6,-20) < (x,y,z) < (69.3,63.6,115), resolved by a 
320 x 300 x 392 grid. In the x-direction (^-direction) , the mesh size is Ax = 0.2 (Ay = 0.2) 
within -20 < x < 20 (-20 < y < 20), and it increases gradually for \x\ > 20 (\y\ > 20), 
reaching about 2 at the boundaries. In the z-direction, the mesh size is Az = 0.2 within 
the range —20 < z < 32, and for z > 32 it increases gradually to 1.39 at the top boundary. 
The horizontal boundaries in the x-direction (the direction of the initial flux tube) are 
assumed to be periodic. The bottom boundary is a perfectly conducting wall and is non- 



penetrating. An outflow boundary condition (see Stone and Norman 1992a) is applied 
for the top boundary and the two horizontal boundaries in the ^-direction. For the initial 
state, we consider a plane-parallel atmosphere in hydrostatic equilibrium, composed of an 
adiabatically stratified polytropic layer representing the top layer of the convection zone, an 
isothermal layer representing the photosphere and the chromosphere, connecting to another 
isothermal layer of 1 MK representing the corona. The vertical profiles of the pressure Po(z), 
density po(z), and temperature T (z) of the hydrostatic field-free atmosphere are shown in 
the bottom panel of Figure 1. A horizontal, twisted magnetic flux tube oriented along the 
x-direction is embedded in the polytropic layer (convection zone) with its axis located at 
(y, z) = (0, —14). The magnetic field of the flux tube is given by 

B = B x (r)± + B e {r)§, (9) 

where 

B x (r) = B ex P (-r 2 /a 2 ), (10) 

B e (r) = qrB x (r). (11) 
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In the above x is the tube axial direction, 9 is the azimuthal direction in the tube cross- 
section, and r denotes the radial distance to the tube axis. The flux tube is uniformly twisted, 
with the constant q denoting the angle of field line rotation about the axis over a unit length 
of the tube. We set B = 5, the radius a = 2, and the rate of twist q = — or 1 (left-handed), 



which makes the tube marginally stable for the kink instability (Linton, Longcope, and 



Fisher |1996 ). The profile of the magnetic pressure p m along the vertical line at (x, y) = (0, 0) 
through the center of the embedded tube is also shown in Figure 1. Inside the flux tube, the 
plasma pressure differs from that of the hydrostatic field- free atmosphere po(z) by: 



Pi(r) 



satisfying the equation: 
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(13) 



which ensures that the pressure gradient in the flux tube balances the Lorentz force. To 
make the flux tube buoyant, we further impose a density change p\ within the flux tube 
relative to po(z) of the field- free hydrostatic atmosphere: 



B 2 



Pi 



-Po[z) 



2po(z) 



[(l + e)exp(-x 2 /A 2 )-e] 



(14) 



where A = 8 and e = 0.2. This makes the middle portion of the flux tube buoyant, with 
Pi/ Po = —0.11 at the middle of the tube axis. The buoyancy declines with \x\ from x = 
following a Gaussian profile with an e-folding length of A, and the two ends are slightly 
anti-buoyant with pi/po = 0.022 at the ends of the tube axis. The plasma (3 (defined as the 
ratio of the plasma pressure over the magnetic pressure) is about 8.3 at the axis of the tube. 



3. Results 
3.1. Overview of the evolution 



The early evolution of the buoyant flux tube is qualitatively similar to those reported 



in many previous investigations (e.g 


Magara and Longcope ||2001 


Fan 


2001 


Magara 


2004 


Manchester et al. 


2004 




Archontis et al. 


2004 




Murray et al. 


2006 


I 


Vlagara 


2006) 



The central portion of the initial flux tube rises buoyantly while the two ends of the flux 
tube sink down towards the bottom owing to the slight negative buoyancy there. The flux 
tube develops into an fi-shape, whose central apex portion rises towards the photosphere. 
When the flux tube enters the stably stratified photosphere, magnetic pressure builds up 
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at the photosphere layer and the magnetic buoyancy instability develops. This causes a 
rapid expansion of the front of the flux tube into the solar atmosphere and initiates the flux 
emergence. 

Figures 2 and 3 show respectively two perspective views of the 3D magnetic field as 
delineated by a set of representative field lines at 4 time instances during the flux emergence 
into the solar atmosphere. We see that at t — 65, the apex of the flux tube has broken through 
the photosphere and the front of the emerging flux tube begins to expand rapidly into the 
atmosphere and corona. With time, a twisted flux rope structure with field lines winding 
about a central axis (as represented by the black line) is formed in the atmosphere (see the 
t = 100 and t = 120 panels). By time t = 120, the central axis of the flux rope structure has 
reached coronal heights, and sigmoid-shaped field lines (see the orange and red field lines), 
some of which containing dips, thread under the central axis with an inverse configuration 
(i.e. pointing from negative to positive polarity region over the neutral line). Note that field 
lines in different images are not the same field lines carrying the same plasma. They are 
simply field lines selectively drawn in each instance to reflect the topological structure of the 
3D field. For example, the black axis in each panel corresponds to a different field line in 
each of the time instances, an important point to be addressed later in §3.2. 

Figures 4(a) (b) show the evolution of the height and rise velocity for 3 different points 
within the central vertical plane of x = 0: (1) the front of the emerging flux tube, (2) 
the Lagrangian position of the original tube axis, and (3) the O-point of the transverse 
magnetic field in the central vertical plane at x — 0, which is the point where B y = on 
the central vertical line of x = and y — 0. Figure 5 shows an example snapshot of the 
central vertical plane of x = at time t = 108, where the arrows show the direction of the 
transverse magnetic field in the plane, overlaid with a color image showing density on a log 
scale. Note, only the part above the photosphere is shown, and the arrow lengths are all 
the same, normalized by the transverse field strength, for the points where the transverse 
field strengths are non-zero. The O-point of the transverse magnetic field is marked by the 
pink '*' point, and it has reached the height of the base of the corona at t — 108. It should 
be emphasized that the O-point position is based on a geometric definition which depends 
on the variation of B y with z along the central vertical line of x = and y — 0. Thus it 
does not track a Lagrangian plasma element, and its height variation with time does not 
in general reflect the true vertical plasma motion but is strongly influenced by the change 
with time of the horizontal magnetic field orientations at points along the central vertical 
line. For these reasons, dz/dt of z(t) for the O-point position shown in Figure 4(a), does not 
in general agree with v z , the vertical velocity of the plasma instantaneously residing at the 
O-point shown in Figure 4(b) (dashed-line). We also note that the black field line shown in 
each of the panels in Figure 2 and 3 is the field line that is drawn from the O-point at that 
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instance, and it corresponds to a different field line in each time. 



From the position and the rise velocity of the front of the flux tube as shown in Figure 
4, one can see that the onset of the magnetic buoyancy instability and flux emergence into 
the atmosphere takes place at about t = 60, at which time there is an exponential increase 



of the vertical speed of the front. Archontis et al. (2004) and Murray et al. (2006) have, 



based on the original work of Newcomb (1961); Thomas and Nye (1975) and Acheson 



(1979), given the following critical condition for the onset of magnetic buoyancy instability 



and flux emergence at the photospheric layer: 



1 + 



k? 



(15) 



where z is height, H p is the local pressure scale height at the photosphere, B is the magnetic 
field strength, 7 is the ratio of the specific heats, (3 is the ratio of the plasma pressure over the 
magnetic pressure, 5 = V — V a d is the superadiabaticity of the atmosphere, which is —0.4 for 
the isothermal stratification of the photosphere, and ku, kj_ and k z are the three components 
of the local perturbation wave vector (normalized by 1/H p ), with k\\ and k± being the two 
horizontal components parallel and perpendicular to the local magnetic field direction, and 
k z being the z component. We find that at t = 60 the build-up of magnetic pressure at the 
photosphere layer has just exceeded the gas pressure, forming a layer of low plasma (3 (see 
top panel of Figure 6), where the strong stabilizing term — {j/2)f35 is reduced to become 
smaller than — H p (d/dz) (log B) (see bottom panel of Figure 6). Using the horizontal widths 
l x and l y of the emerging magnetic region at z = 2 to estimate k\\ ~ 2i\ jly and k± ~ 2%/l x , 
and k z ~ 2n / tube radius, we find that the above critical condition is met at z ~ 2. 

As the front of the emerging flux tube begins to expand into the atmosphere at about 
t = 60, the rise speed of the original tube axis (solid curve in Figure 4(b)) and the rise speed 
measured at the location of the O-point (dashed curve in Figure 4(b)) both begin to decline, 
and their positions remain below z = 5 during the period from t = 65 to t = 86 (see Figure 
4(a)). However, soon after t = 86 there is an upward jump of the O-point position, and then 
its height continues to increase sharply. With the O-point moving into coronal heights, a 
flux rope with dipped, sigmoid-shaped core field (see the t = 120 panels in Figures 2 and 3) 
develops in the corona. 



3.2. Vortical motion of the two polarities and flux rope formation 



To understand the apparent rise of the O-point into coronal heights and the formation 
of the coronal flux rope, we have tracked the evolution of a set of field lines which intersect 
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the center vertical line of the domain (at x — and y — 0). The way we do this is by tracking 
the motion of a few Lagrangian points situated along the center vertical line of symmetry at 
x = and y — 0, and trace out field lines through these points at successive time instances 
from t = 65 to t = 120. Figures 7, 8, and 9 show 3D evolution of these tracked field lines 
as viewed from 3 perspectives. The black field line here corresponds to the original tube 
axis, and all the other field lines have their mid points (at x = 0, y = 0) above the mid 
point of the black field line (reddish field lines have mid points above bluish field lines). 
Note that the apex of the red field line does not correspond to the emerging front, but is 
significantly below it. Thus its position and rise speed do not correspond to the dash-dotted 
lines for the front of the emerging tube shown in Figure 4. In this sequence of images, a 
field line of a particular color corresponds to the same field line carrying the same plasma. 
A corresponding animated GIF movie showing the evolution of the tracked field lines is also 
available from the electronic edition. 

At the beginning of the flux emergence (see the t = 65 image in Figure 8), the field 
lines all show a normal configuration, i.e. they arch over the polarity inversion line from 
the positive to the negative polarity region. The B y components of these field lines at the 
midpoints are all positive. Thus at this time, the O-point, which is the point of zero B y 
on the central vertical line, is located below the black field line (the original tube axis). 
With time the two polarity regions on the photosphere begin to show a significant shearing 



motion along the neutral line (e.g. 


Manchester 


2001 


Manchester et al. 


2004 


Fan 


2001 


Magara and Longcope 


2003; Magara 


2004, 2006 


), and the two polarity regions begin to 



separate. As a result the lower-lying black (original axial field line) and dark blue field lines 
are lengthened and begin to rotate (clockwise) toward an inverse configuration, with their 
By changing from positive to negative at the midpoints (see Figure 8). By time t = 90 (see 
Figure 8), the black and dark blue field lines have rotated to show a slight negative B y at 
their mid points, while the upper yellow and red field lines are still having positive B y at 
their midpoints. This change of orientations of the lower (black and dark blue) field lines 
from positive B y to negative B y causes an abrupt upward shift of the zero point of B y along 
the center vertical line of x = and y = 0, i.e. an upward shift of the O-point seen in Figure 
4(a) (soon after t = 85). As a result the O-point position shifts from below to above the 
original tube axis (the black field line in Figures 7,8 and 9). This apparent crossing of the 
O-point position from below to above the original tube axis (as seen between t = 85 and 
t = 90 in Figure 4a) is not a true plasma rise motion, but is simply due to the change of 
horizontal field line orientations in the atmosphere: the original tube axis (the black field 
line tracked in Figure 8) rotates from a forward into an inverse configuration, with B y at its 
mid point changing from positive to negative, and as a result the zero point of B y on the 
center vertical line of x = and y = 0, namely the O-point, propagates upward in position 
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to a different field line above, which is now aligned with the x direction. 

By time of about t = 90, a significant counter-clockwise vortical or rotational motion 
has set in within each of the two polarity flux concentrations on the photosphere, reminiscent 



of sunspot rotations (e.g. Brown et al. 2003). The vortical motions can be seen from the 



footpoint motions of the field lines in Figures 7,8, and 9 (and in the movie), and more clearly 
in Figure 10 which shows snapshots of the z component of the vorticity uo z on the photosphere 
overlaid with contours of B z . We find high intensity of positive u z (corresponding to counter- 
clockwise vortical motion) centered on the peaks of the vertical magnetic flux concentration 
of both polarities. Figure 11 shows the evolution of the mean vertical vorticity < uj z > 
averaged over the area of each polarity flux concentration where B z is above 75% of the peak 
B z value. We find that the counter-clockwise vortical motion centered on each polarity flux 
concentration, as represented by < u z >, rises to its peak rate at t ~ 93 and then persists at 
a relatively steady rate throughout the subsequent evolution. The footpoint vortical motions 
further twist up the emerged field lines, causing continued clock-wise horizontal rotation of 
the blue and the higher yellow and red field lines in the atmosphere as viewed from the 
top (see the t = 100, t = 110, t = 118, and t = 120 panels in Figure 8 and the movie). 
With time, the yellow and then red field lines are also progressively rotated into an inverse 
configuration, changing the sign of B y from positive to negative at their mid points. This 
causes the O-point position (the zero B y position on the center line of x = 0, y = 0) to 
continue to propagate upward as seen in Figure 4(a) after about t — 90. Note again that 
this upward shift of the O-point position from t = 86 to about t = 102 is not due to a true 
vertical rise of the plasma as can be seen from the near zero v z during this period shown 
in Figure 4(b), but is simply a propagation effect produced by the change of horizontal 
orientations of the emerged field lines. 

Starting from about t = 102, the plasma rise velocity v z measured at the O-point 
position also shows a significant upward acceleration (see Figure 4(b)), which means that 
the rapid ascent of the O-point after about t = 102 (see Figure 4(a)) is a combination of the 
upward propagation produced by the change of the horizontal orientations of the field lines 
and an actual rise motion of the field lines. At about t = 118 (see Figure 8), the O-point 
position has reached z = 40.2, which is at this instant identified with the mid-point of the 
red field line. The red field line has rotated to become aligned with the x-direction with zero 
By at the middle and is at this instant identified as the axial field line of the coronal flux 
rope. The plasma rise speed v z measured at the O-point is about 2.6 (or 16.9 km/s), which 
continues to increase (see Figure 4(b)). The field lines continue to rotate horizontally. As a 
result the field line identified as the coronal axial field line is continually changing. By time 
t = 120, the red field line has further rotated to having a negative B y at its center, and the 
O-point has now shifted to a field line (not shown in Figures 7, 8, and 9) above the red field 
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line. 

We also find that in conjunction with the horizontal rotation of the black and blue 
field lines in the atmosphere due to shearing and twisting at their footpoints, they develop 
central dips, into which plasma drains, and the two side "wings" of each of these field lines 
rise rapidly into the corona (see the t = 100, t = 110, t = 118, t = 120 images in Figures 7 
and 9). 

The above examination of the tracked field lines indicate that the formation and rise 
of the coronal flux rope is not a result of the bodily emergence of the subsurface flux tube 
as a whole. The bottom concave portions of the helical field lines in the original emerging 
flux tube breaking through the surface remain largely trapped at and below the photosphere 
layer due to the weight of the heavy plasma, and the center of the original tube axis ceases 
to rise after reaching about 2-3 pressure scale heights above the photosphere. However, the 
shearing and vortical motions of the two polarity flux concentrations at the photosphere 
continue to transport twist and magnetic energy into the emerged field in the atmosphere 



e.g Magara and Longcope 2003 Manchester et al. 2004 Magara 2006). Our analysis in 



this paper shows that due to this vortical motion at the footpoints, the emerged field lines 
near and just above the original axial field line (including the original axial field line itself) 
rotate horizontally in the atmosphere, changing B y at their mid points from positive to 
negative (i.e. from an original normal configuration to an inverse configuration with respect 
to the photosphere neutral line). This initiates an upward propagation of the O-point in the 
emerged field in the central x = cross-section simply because of the change of horizontal 
orientation of the inner field lines, and not because of the rise of these field lines. Later, 
the plasma rise speed measured at the O-point also shows a sharp increase, contributing to 
the rapid ascent of the O-point into coronal heights. As a result, a coronal flux rope with 
an axis identified with the field line instantaneously aligned with the x-direction at its mid 
point (i.e. going through the O-point), and with sigmoid shaped field lines threading under 
the axis with an inverse configuration develops, and it is rising with increasing velocity as 
the vortical motions at the footpoints of the flux rope continue to transport magnetic energy 
and twist into the corona. Note that the field line being identified as the axial field line of 
the coronal flux rope is continually changing, due to the continued horizontal rotation of the 
field lines in the atmosphere as their footpoints are being twisted. 

To quantify the transport of twist into the solar atmosphere due to both flux emergence 
as well as the horizontal vortical motions at the photosphere, we have computed the evolution 
of the relative magnetic helicity in the atmosphere above the photosphere and the helicity 
flux through the photospheric boundary. The relative magnetic helicity in the atmosphere 



- 12 - 



above the photosphere z = is defined as ( Berger and Field~||1984 Finn and Antonsen"|T9 85): 



Hrr 



(A + A v ) • (B - P) dV 



(16) 



z>0 



where B is the magnetic field in the unbounded half space above z = 0, A is the vector 
potential for B, P is the reference potential field having the same normal flux distribution 
as B on the z = boundary, and A p is the vector potential for P. The relative magnetic 



helicity H m is invariant with respect to the gauges for A and A p . Following Devore (2000), 
we use the following specific A and A p for computing H m : 



(17) 



where, 



A(x,y,z) = A p (x,y,0) - z x / B(x,y,z')dz' 

Jo 

POO 

A p (x, y, z) = V x z / 0(x, y, z) dz f , 

J z 

B z (x' : y',0) 



(?,y,z) 



2tt 



[(x - x'Y + (y - y') 2 + z 2 } 1 / 



Yn dx d V 



(19) 



If we use the specific A and A p given in equations (17) through (19), it can be shown that 



the expression (16) for the relative helicity H m reduces to simply 



A ■ B dV. 



(20) 



z>0 



Note, the derivation of the relative magnetic helicity assumed unbounded half space above 



the photosphere (Berger and Field 1984 Devore 2000). We carry out the above integral 



within the simulation domain above z = assuming that the field outside of the domain is 
zero. Thus, our helicity integration only gives an accurate result before the emerged field has 



expanded to reach the top and side boundaries. It has also been shown (Berger and Field 



1984 Devore 2000 ) that the rate of change of the relative helicity H m can alternatively be 
computed by evaluating the helicity flux through the z = boundary plane: 



dH n 
~dt 



-2 J J [(A p ■ v)B - (A p ■ B)v] • z dxdy 



(21) 



where A p is given by equations ( 18 )-( 19 ). The first term on the right-hand-side of equation 



(21) corresponds to the injection of helicity by horizontal motion on the boundary, and the 



second term corresponds to injection by direct emergence or submergence of twisted flux 
across the boundary. Note again, we carry out the above integration on the z = plane 
within the domain, assuming zero flux through the z = plane outside of the domain. The 
resulting evolution of dH m /dt and H m for the emerged field above z = is shown in Figure 



13 



12. In the lower panel of Figure 12, both H m computed from equation (16) (solid line) and 
that obtained by accumulation of helicity flux with time H m = J (dH m /dt) dt (dash-dotted 
line) are shown. The good agreement of the two indicates good helicity conservation from 
our simulation. The deviation in the end is due to the fact that twisted magnetic flux is 
exiting the top and side boundaries and is therefore not accounted for in the integration of 



equation 16 while the time integration of dH m /dt continues to give accurately the amount of 
helicity that has been transported into the atmosphere through the boundary at z — 0. The 
top panel of Figure 12 shows that the helicity flux due to flux emergence is dominating for 
a brief period at the beginning (from about t — 55 to t — 80), and then after about t = 80, 
horizontal shearing and rotational motions at the footpoints of the emerged field become the 
main source of helicity flux. By time of about t = 90, the rotational motions in the two 
polarity flux concentrations have been established (see Figure 11) and they provide a steady 
and dominant helicity flux into the atmosphere in the subsequent evolution (top panel of 
Figure 12). By time t = 120, the helicity transported into the atmosphere has reached about 
— $^ ube , corresponding to one full twist (2n rotation) of the flux $ tu b e in the original flux 
tube. 

The rotational motions centered on the two polarity flux concentrations are a manifes- 
tation of non-linear torsional Alfven waves propagating along the flux tube, transporting 
twist from the tube's interior portion towards its expanded coronal portion. Consider the 
quantity a = (VxB)-B/B 2 , which gives a measure of the local rate of twist of the magnetic 
field. For a force-free magnetic field: 

V x B = aB, (22) 

from which one can show that 

B ■ Va = 0, (23) 

i.e. a is constant along each field line for a field with zero Lorentz force. If there is a 
gradient of a along a field line, then the thin flux bundle (or tube) centered on the field line 



will experience a net axial torque which drives torsional motion of the flux tube (Longcope 



and Klapper 1997). The net axial torque exerted on a thin flux tube segment (1,1 + dl), 



where I denotes the arc length along the thin tube, is given by Longcope and Welsch (2000 
eq. (20) in that paper): 

Tl = ±a A B?^dl (24) 

' 16 1 dl y ' 

where T\ represents a torque along the axial direction 1, which will drive the tube segment 
to spin about its axis, a denotes the radius of the tube, Bi denotes the axial field strength. 



Equation (24) states that the net axial torque T\ for a thin flux tube segment is directly 



proportional to the gradient of a along the tube. Note that in Longcope and Welsch (2000), 
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the quantity q, denned as the angle of field line rotation about the central axis per unit length 
along the tube, is used in place of a to describe the twist. For a thin flux tube, the relation 
between a and q is a = 2q (Longcope and Klapper 1997). Figure 13 shows the variation of 
a along three field lines: the black field line shown in Figures 7, 8, and 9, corresponding to 
the original tube axis, and its two neighboring blue field lines shown in Figures 7, 8, and 9, at 
t — 118. We have plotted the a values along these three field lines (using the same colors for 
the data points as those of the corresponding field lines in Figures 7, 8, and 9) as a function 
of depth, from their left ends to the left apices in the atmosphere, at t — 118. We can see 
that for all three neighboring field lines, the magnitude of a decreases with height from the 
interior into the atmosphere. The a is nearly constant along the atmosphere/corona portions 
of the field lines, indicating that they are near a force-free state, and the coronal value is 
about a factor of 10 less than the peak magnitude of a in the interior part of the flux tube. 
The low a magnitude in the corona is due to the extreme expansion and stretching of the 
magnetic field lines resulting from flux emergence (Longcope and Welsch |2000| >. The lowered 
magnitude of a in the atmosphere establishes a gradient of a along the flux tube from the 
interior to the atmosphere, and produces a torque that drives the observed torsional motion 
of the flux tube intersecting the photosphere. The torsional motion will continue until the 
gradient of a along the tube is removed, and the coronal a equilibrates with the interior a 
(|Longcope and Welsch |[2000|. 



3.3. Sigmoid fields and current sheet: 



The final 3D magnetic field in the atmosphere and corona obtained from the simulation 
show a coronal flux rope structure whose properties may qualitatively explain some basic 



observed features of sigmoid active regions. Similar to what was found in (Manchester et al. 



2004 


Magara 


2004, 


2006) 



panel of Figure 14 shows the current density J in a horizontal plane at height z = 10 in the 
chromosphere. It shows a horizontal cross-section of the current sheet. The current sheet 
has an over-all inverse-S shape (except the central zigzag). Field lines traced through a few 
points along the current concentration on z = 10 are also plotted (for different perspectives) 
in the two lower panels of Figure 14. We find that all the field lines going through the current 
sheet show an inverse-S shape when viewed from the top (as illustrated in the lower right 



panel of Figure 14), similar to what has been found in the simulations of Magara (2004 



2006) and consistent with the theoretical prediction of Titov and Demoulin (1999) and Low 



and Berger (2003). One can argue that the heating associated with the current sheet may 
cause these inverse-S shaped field lines to preferentially brighten up in soft X-ray, giving rise 
to the observed X-ray sigmoid loops in an active region. It can also be seen that some of 
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the field lines undergo sharp dips at the center portion of the current sheet (see the red field 
lines in the lower left panel of Figure 14). Reconnections in the current sheet are disengaging 
the heavy mass trapped in the dips, allowing for the flux rope to rise further into the corona 
with increasing speed. 

We find that the inverse-S shaped current sheet begins to develop at the time of about 
t = 102, being initially confined to a few pressure scale-heights above the photosphere. With 
time the current sheet extends upward, and by time t = 120, it has extended to the base 
of the corona. Growing reconnections in the current sheet due to the numerical diffusion 
inherent in the code are expected to aid the upward acceleration of the center of the flux 
rope seen after about t = 102 (Figure 4b). However, the ultimate driver of the rise of the 
flux rope and the development of the current sheet is the continued vortical motions at the 
footpoints of the emerged field lines, which transport twist and magnetic energy into the 
atmosphere, causing an excess of the upward Lorentz force that drives the acceleration. 



Discussions 



We have presented a 3D MHD simulation of the emergence of a twisted f2-shaped flux 
tube from the top layer of the solar interior into the solar atmosphere and the corona. 



Compared to some of the previous simulations (e.g. Fan 2001 Archontis et al. 2004 



Manchester et al. 2004), the subsurface rising flux tube is more sharply arched such that 



the emerging segment first breaking through the photosphere contains less than one full 
wind of twist. In this simulation, we find that after an initial stage of flux emergence, during 
which the two polarity flux concentrations on the photosphere undergo a shearing motion 
and become separated, a prominent rotational/vortical motion develops within each polarity, 
reminiscent of sunspot rotations. This rotational motion persists throughout the subsequent 
evolution and steadily transports twist or magnetic helicity into the atmosphere. 



Similar to what has been found in many previous investigations (e.g. Fan ||2001 Archon 



tis et al. 2004 Manchester et al. 2004 Magara 2006), the twisted flux tube is unable to 



rise bodily into the corona as a whole. While the upper parts of the winding field lines of the 
twisted flux tube expand into the atmosphere and the corona, the center of the original axial 
field line of the tube ceases to rise at about 2-3 pressure scale heights above the photosphere 
and the U-shaped parts of the helical field lines largely remain anchored at and below the 
photospheric layer. However, we find that the initial shearing and the subsequent rotational 
motions of the two polarity flux concentrations twist up the inner field lines (near and in- 
clude the original tube axis) of the emerged field, causing them to rotate in the atmosphere 
from their initial normal configuration (having a positive B y at the mid point, or arching 
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over the polarity inversion line from the positive to the negative polarity), into an inverse 
configuration (having a negative B y at the mid point, directed from the negative polarity 
to the positive polarity over the neutral line). This change of the horizontal orientations of 
the emerged field lines causes an upward propagation of the O-point (the zero-crossing point 
of B y along the center vertical line of x = and y = 0) in the emerged field in the central 
vertical cross-section (x = 0). With the continued twisting of the emerged field lines by 
the vortical motions at the footpoints, which transports twist and magnetic energy into the 
corona, the actual plasma motion measured at the O-point also starts to accelerate upwards. 
Thus, the rapid ascent of the O-point is a combination of the upward propagation produced 
by the change of the horizontal orientations of the field lines and the actual rise of the field 
lines. As a result, a coronal flux rope with an axial field line identified with the field line in- 
stantaneously aligned with the x-direction at its mid point (i.e. going through the O-point), 
and with sigmoid shaped field lines threading under the axis with an inverse configuration, 
is formed and rises with increasing speed as the vortical motions at the footpoints continue. 
Some of the sigmoid shaped field lines threading under the axis have also develop sharp dips 
as they are being twisted at the footpoints, with plasma draining into the dips and with the 
two side lobes rising rapidly into the corona. A current sheet forms in the lower atmosphere 
with field lines going through it all showing sigmoid shapes. Heating in the current sheet 
may therefore cause these field lines to brighten up in X-ray as the observed sigmoid loops 
in active regions. Field lines going through the center portion of the current sheet undergo 
sharp dips. Reconnections in the current sheet allow for removal of the heavy plasma trapped 
in the dips and hence enhance the upward acceleration of the coronal flux rope. 

The formation of a coronal flux rope with sigmoid shaped core fields and with a current 
sheet in the lower atmosphere as a result of the emergence of a twisted flux tube from the 



interior into the solar atmosphere has been found to be a very general result (e.g. Manch 
ester et al. 1120041 iMagara 1120041 120061 lArchontis and Toeroek 1 120081). The onset of a 



rotational/vortical motion centered on the two polarity flux concentrations soon after they 



become separated has also been found in Magara (2006). Our work explicitly demonstrates 



how the vortical motions of the footpoints twist up the inner emerged field lines and rotate 
them into an inverse configuration (Figures 7, 8, and 9), leading to the development and 
upward ascent of the coronal flux rope structure. We also show that the rotational motions 
of the two polarity flux concentrations on the photosphere are a manifestation of nonlinear 
torsional Alfven waves propagating along the flux tube, consistent with what has been pre- 



dicted by an idealized analytical model of Longcope and Welsch (2000). Due to the rapid 



stretching of the emerged magnetic field in the atmosphere and corona, the magnitude of 
a = J ■ B/5 2 along the coronal field lines, which is a measure of the rate of twist per unit 
length, drastically decreases. As a result, a gradient of the rate of twist is established along 
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the flux tube from the interior into the atmosphere with the interior portion of the flux 
tube having a much higher rate of twist. This gradient in the rate of twist drives torsional 
waves along the flux tube, transporting twist from the interior highly twisted portion into 



the expanded coronal portion (Longcope and Welsch 2000). Thus the rotational motion 



will continue until the coronal a equilibrates with the interior a along the field lines. The 
time scale for establishing the equilibrium is on the order of the Alfven transient time along 
the interior flux tube, which means that the rotational motion can persists for a few days 
after the initial emergence. Depending on how great the stretching of the emerged field 
lines is compared to the total length of the interior flux tube, it is conceivable that the final 
equilibrated a value is of a significantly smaller magnitude compared to the initial a of the 
subsurface flux tube rising to the photosphere. This may provide an explanation for the 
significantly lower a values measured for the majority of the solar active regions compared 
to the twist rate needed for a buoyant flux tube to rise cohesively in the solar interior. 

Sunspot rotations have been observed in many events preceding X-ray sigmoid bright- 
ening and the onset of eruptive flares (e.g Brown et al.||2003 ; Zhang et al.||2008 ). Our results 



suggest that these rotational motions are due to non-linear torsional Alfven waves naturally 
occurring during the emergence of a twisted flux tube, and is an important means whereby 
twist is transported from the interior into the corona, driving the development of a coronal 
flux rope as a precursor structure for solar eruptions. The horizontal vortical motion and 
its subsurface extension corresponding to the torsional Alfven waves may be detectable by 
local helioseismology techniques. 
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Fig. 1. — The initial setup of the simulation. The top panel shows the simulation domain, 
where the green plane denotes the photosphere and the configuration of the initial flux tube 
is illustrated with 3 traced field lines in the flux tube. The bottom panel shows the vertical 
profiles of plasma pressure, density and temperature of the field-free static atmosphere, and 
also the profile of the magnetic pressure through the embedded flux tube along the center 
line at (x, y) = (0, 0). 
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t=65. t=80. 




Fig. 2. — Snapshots of the 3D magnetic field during the flux emergence into the atmosphere. 
Note that field lines in the different snapshots are not the same field lines carrying the 
same plasma - they are simply field lines selectively drawn in each instance to reflect the 
topological structure of the 3D field. The black field line shown in each of the panels is the 
field line that is drawn from the O-point of the transverse field in the central cross-section 
(at x = 0), and it corresponds to a different field line in each time. 




Fig. 3. — Same as Figure 2 but viewed from a different perspective. 
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Fig. 4. — The evolution of the height (a) and vertical velocity (b) for 3 different points within 
the central vertical plane of x = 0: (1) the front of the emerging flux tube (dash-dotted line), 
(2) the Lagrangian position of the original tube axis (solid line), and (3) the O-point of the 
transverse magnetic field in the central cross-section (at x — 0), which is the point where 
By = on the center vertical line of x = and y — 0. The mid-points of the black field lines 
shown in Figure 2 and 3 correspond to the O-point positions at those instances. 
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t=108.0 




Fig. 5. — A snapshot of the central vertical plane of x — at t — 108, with arrows showing the 
direction of the transverse magnetic field (lengths of the arrows are all the same, normalized 
by the transverse magnetic field strength for all points where the transverse field strength is 
non-zero), overlaid with a color image showing the density in log scale. Only the part above 
the photosphere is shown. The O-point of the transverse magnetic field is marked by the 
pink '*' point 
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Fig. 6. — Vertical profiles of several quantities near the photosphere along the center line 
of (x,y) = (0,0), at t = 60. It shows the conditions of the emerging field entering the 
photosphere when it becomes unstable to the onset of the magnetic buoyancy instability 
(see discussion in the text) 
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Fig. 7. — 3D evolution of a set of tracked field lines as they are being twisted up by the shear 
and rotation motions at their footpoints on the photosphere. In this sequence of images, a 
field line of a particular color corresponds to the same field line carrying the same plasma. 
The black field line corresponds to the original tube axis, and all the other field lines have 
their mid points (at x = 0, y — 0) above the mid point of the black field line (reddish field 
lines have mid points above bluish field lines). An animated GIF movie is also available in 
the electronic version. 




Fig. 8. — Same as Figure 7, but is viewed from the top 



- 28 - 



t=S5,0 t=90.0 




Fig. 9. — Same as Figure 7, but a different side view perspective. 
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t= 90.01 




Fig. 10. — Snapshots of the z component of the vorticity oo z on the photosphere overlaid 
with contours of B z with solid (dashed) contours representing positive (negative) B z . The 
images show counter-clockwise vortical motion centered on the peaks of the vertical flux 
concentration of both polarities. 
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Fig. 11. — The evolution of the mean vertical vorticity < uo z > averaged over the area of 
each polarity flux concentration where B z is above 75% of the peak B z value. 
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time 

Fig. 12. — Top panel shows the helicity injection into the atmosphere dH m /dt due to hor- 
izontal motions on the photosphere (dash-dotted line), vertical emergence of flux (dashed 
line), and the sum of the two (solid line). The bottom panel shows the evolution of the 
relative magnetic helicity H m of the emerged field computed from equation (6) (solid line) 
and through time integration of the total dH m /dt (dash-dotted line). 
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Fig. 13. — The variation of a = (V X B) ■ ~B/B 2 as a function of z along three field lines: 
the black field line shown in Figures 7, 8, and 9, which is the original tube axis, and its two 
neighboring blue field lines shown in Figures 7, 8, and 9, at time t = 118. The a values are 
plotted along these three field lines (using the same colors for the data points as those of the 
corresponding field lines in Figures 7, 8, and 9) as a function of depth, from their left ends 
to the left apices in the atmosphere. 
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Fig. 14. — The top panel shows the current density J = V x B in a horizontal plane at 
height z = 10 in the chromosphere. The two lower panels show two perspective views of a 
set of 3D field lines traced from a few points along the current concentration shown in the 
top panel. 



